GDP per capita

This is GDP per capita, PPP (constant 2021 international $). Retrieved from the World Bank (https://data.worldbank.org/indicator/NY.GDP.PCAP.PP.KD) for the year 2019.

gdp_data <- read_csv("data/gdp_pc_ppp_2021.csv", skip = 4, show_col_types = FALSE)

head(gdp_data)

## # A tibble: 6 × 69
##   `Country Name`  `Country Code` `Indicator Name` `Indicator Code` `1960` `1961`
##   <chr>           <chr>          <chr>            <chr>            <lgl>  <lgl> 
## 1 Aruba           ABW            GDP per capita,… NY.GDP.PCAP.PP.… NA     NA    
## 2 Africa Eastern… AFE            GDP per capita,… NY.GDP.PCAP.PP.… NA     NA    
## 3 Afghanistan     AFG            GDP per capita,… NY.GDP.PCAP.PP.… NA     NA    
## 4 Africa Western… AFW            GDP per capita,… NY.GDP.PCAP.PP.… NA     NA    
## 5 Angola          AGO            GDP per capita,… NY.GDP.PCAP.PP.… NA     NA    
## 6 Albania         ALB            GDP per capita,… NY.GDP.PCAP.PP.… NA     NA    
## # ℹ 63 more variables: `1962` <lgl>, `1963` <lgl>, `1964` <lgl>, `1965` <lgl>,
## #   `1966` <lgl>, `1967` <lgl>, `1968` <lgl>, `1969` <lgl>, `1970` <lgl>,
## #   `1971` <lgl>, `1972` <lgl>, `1973` <lgl>, `1974` <lgl>, `1975` <lgl>,
## #   `1976` <lgl>, `1977` <lgl>, `1978` <lgl>, `1979` <lgl>, `1980` <lgl>,
## #   `1981` <lgl>, `1982` <lgl>, `1983` <lgl>, `1984` <lgl>, `1985` <lgl>,
## #   `1986` <lgl>, `1987` <lgl>, `1988` <lgl>, `1989` <lgl>, `1990` <dbl>,
## #   `1991` <dbl>, `1992` <dbl>, `1993` <dbl>, `1994` <dbl>, `1995` <dbl>, …

gdp_data <- gdp_data |> 
  select(`Country Name`, `Country Code`, `2019`)


gdp_data <- gdp_data |> 
  rename(country_name = `Country Name`,
         country_code = `Country Code`,
         gdp_pc_ppp = `2019`)

gdp_data <- gdp_data |> 
inner_join(select(codelist, iso3c), by = c("country_code" = "iso3c"))

str(gdp_data)

## tibble [28 × 3] (S3: tbl_df/tbl/data.frame)
##  $ country_name: chr [1:28] "Austria" "Belgium" "Bulgaria" "Cyprus" ...
##  $ country_code: chr [1:28] "AUT" "BEL" "BGR" "CYP" ...
##  $ gdp_pc_ppp  : num [1:28] 64630 60452 27673 46157 47720 ...

We now have 28 observations including the UK.

LGBT rights

Found this LGBT rights index on Our World in Data (https://ourworldindata.org/grapher/lgbt-rights-index) which captures whether LGBT+ people enjoy the same rights as cisgender people combining information on 18 different policies. It includes the legal status of same-sex marriage.

lgbt_rights_index <- read.csv("data/lgbt-rights-index.csv")

lgbt_rights_index <- lgbt_rights_index |> 
  filter(Year == 2019) |> 
  inner_join(select(codelist, iso3c), by = c("Code" = "iso3c")) |> 
  select(-Year) # all observations are from 2019

str(lgbt_rights_index)

## 'data.frame':    28 obs. of  3 variables:
##  $ Entity            : chr  "Austria" "Belgium" "Bulgaria" "Croatia" ...
##  $ Code              : chr  "AUT" "BEL" "BGR" "HRV" ...
##  $ LGBT..Policy.Index: num  8.92 10 4.94 6.96 3.95 ...

Column names can be renamed:

lgbt_rights_index <- lgbt_rights_index |> 
  rename(country_name = Entity,
         country_code = Code,
         lgbt_policy_index = LGBT..Policy.Index) 

Gender inequality index

Gender Development and Gender Inequality indexes, developed by the United Nations Development. Retrieved from https://hdr.undp.org/data-center/documentation-and-downloads for the year 2019.

gender_index <- read_xlsx("data/UNDP_gender_indexes.xlsx")

gender_index <- gender_index |>
  select(-dimension, -note, -year) |> # empty/useless columns
  inner_join(select(codelist, iso3c), by = c("countryIsoCode" = "iso3c"))

gender_index |> 
  distinct(country) |> 
  nrow()

## [1] 28

str(gender_index)

## tibble [560 × 7] (S3: tbl_df/tbl/data.frame)
##  $ countryIsoCode: chr [1:560] "AUT" "AUT" "AUT" "AUT" ...
##  $ country       : chr [1:560] "Austria" "Austria" "Austria" "Austria" ...
##  $ indexCode     : chr [1:560] "GII" "GDI" "GDI" "GDI" ...
##  $ index         : chr [1:560] "Gender Inequality Index" "Gender Development Index" "Gender Development Index" "Gender Development Index" ...
##  $ indicatorCode : chr [1:560] "abr" "eys_f" "eys_m" "gdi" ...
##  $ indicator     : chr [1:560] "Adolescent Birth Rate (births per 1,000 women ages 15-19)" "Expected Years of Schooling, female (years)" "Expected Years of Schooling, male (years)" "Gender Development Index (value)" ...
##  $ value         : num [1:560] 5.499 16.451 15.725 0.969 0.053 ...

We checked that we have indeed 28 distinct countries in the dataset.

This dataset collects many indicators apart from the indexes values. We are selecting only the Gender Development and Gender Inequality indexes (GDI and GII):

gender_index <- gender_index |> 
  filter(indicatorCode %in% c("gdi", "gii")) |> 
  select(-c(indexCode, indicatorCode, indicator)) |> # removing redundant columns
  pivot_wider(names_from = index, values_from = value)

colnames(gender_index) <- janitor::make_clean_names(colnames(gender_index)) # cleaning spaces and upper cases

And then exploring them to see which one is a better fit for modeling and predictions:

summary(gender_index$gender_inequality_index)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0160  0.0555  0.0890  0.1067  0.1335  0.2390

summary(gender_index$gender_development_index)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.9560  0.9725  0.9855  0.9856  0.9918  1.0230

sd(gender_index$gender_inequality_index)

## [1] 0.06664423

sd(gender_index$gender_development_index)

## [1] 0.0167627

We observe that the second one (Gender Development Index) has a very small range, indicating that most countries have nearly identical scores. This is supported by the standard deviation, which shows that the GII is more spread out than the GDI. Being almost constant, GDI won’t add much value to the analysis. This is likely due to the fact that these indexes are created by the United Nations Development Programme for all countries, so it may not capture the finer differences between EU countries.

cor(gender_index$gender_inequality_index, 
    gender_index$gender_development_index, 
    use = "complete.obs")

## [1] 0.1750589

Surprisingly, the two indexes have a very weak positive correlation. While they are not inverse, we interpreted them as measuring opposite things, and had expected a negative correlation. Maybe the correlation coefficient is not useful here because of the low variation of both variables (specially GDI), so ultimately we chose to keep GII and discard GDI.

gender_index <- gender_index |> 
  select(-gender_development_index)

Economist’s Democracy Index

democracy_index <- read_xlsx("data/EIU_democracy_index.xlsx", sheet = 4)

# the ISO codes were lowercase which impedes the join
democracy_index$geo <- toupper(democracy_index$geo)

# filter for 2019 and EU28 countries
democracy_index <- democracy_index |> 
  filter(time == 2019) |> 
  inner_join((select(codelist, iso3c)), by = c("geo" = "iso3c"))

# clean var names
names(democracy_index) <- names(democracy_index) %>%
  janitor::make_clean_names() %>%
  gsub("_eiu$", "", .)

democracy_index <- democracy_index |> 
  rename(country_code = geo,
         country_name = name,
         year = time)

str(democracy_index)

## tibble [28 × 10] (S3: tbl_df/tbl/data.frame)
##  $ country_code                 : chr [1:28] "AUT" "BEL" "BGR" "HRV" ...
##  $ country_name                 : chr [1:28] "Austria" "Belgium" "Bulgaria" "Croatia" ...
##  $ year                         : num [1:28] 2019 2019 2019 2019 2019 ...
##  $ democracy_index              : num [1:28] 82.9 76.4 70.3 65.7 75.9 76.9 92.2 79 92.5 81.2 ...
##  $ electoral_pluralism_index    : num [1:28] 95.8 95.8 91.7 91.7 91.7 95.8 100 95.8 100 95.8 ...
##  $ government_index             : num [1:28] 78.6 82.1 64.3 60.7 64.3 67.9 92.9 78.6 89.3 78.6 ...
##  $ political_participation_index: num [1:28] 83.3 50 72.2 55.6 66.7 66.7 83.3 66.7 88.9 77.8 ...
##  $ political_culture_index      : num [1:28] 68.8 68.8 43.8 50 68.8 68.8 93.8 68.8 87.5 68.8 ...
##  $ civil_liberties_index        : num [1:28] 88.2 85.3 79.4 70.6 88.2 85.3 91.2 85.3 97.1 85.3 ...
##  $ change_in_democracy_index    : num [1:28] 0 -1.4 0 0 0 ...

In this case, we are only keeping the overall index value.

democracy_index <- democracy_index |> 
  select(country_name, country_code, democracy_index)

Joining all country-level data together

country_level_data <- codelist |> 
  select(iso3c, iso2c) |> 
  left_join(gdp_data, by = c("iso3c" = "country_code")) |> 
  # left_join(rural_data, by = c("iso3c" = "country_code")) |> rural population data will probably be discarded
  left_join(gender_index, by = c("iso3c" = "country_iso_code")) |> 
  left_join(lgbt_rights_index, by = c("iso3c" = "country_code")) |> 
  left_join(democracy_index, by = c("iso3c" = "country_code")) |> 
  select(-contains("country_")) # removing all duplicated country_name columns that were joined from each data frame
str(country_level_data)

## tibble [28 × 7] (S3: tbl_df/tbl/data.frame)
##  $ iso3c                  : chr [1:28] "AUT" "BEL" "BGR" "HRV" ...
##  $ iso2c                  : chr [1:28] "AT" "BE" "BG" "HR" ...
##  $ gdp_pc_ppp             : num [1:28] 64630 60452 27673 35094 46157 ...
##  $ country                : chr [1:28] "Austria" "Belgium" "Bulgaria" "Croatia" ...
##  $ gender_inequality_index: num [1:28] 0.053 0.048 0.205 0.119 0.235 0.124 0.016 0.092 0.031 0.086 ...
##  $ lgbt_policy_index      : num [1:28] 8.92 10 4.94 6.96 3.95 ...
##  $ democracy_index        : num [1:28] 82.9 76.4 70.3 65.7 75.9 76.9 92.2 79 92.5 81.2 ...

Extract data labels

Exploiting the fact that the .dta files has attributes (labels) for all its columns.

If we search for attr(colname, "label") you get back the original long name of the variable

If we instead search for attr(colname, "labels") you get back all the possible encoding levels of the variable (ex.1,2,3,4 etc.) and by doing names(attr(colname, "labels")) you get back the actual meaning of those numbers (ex. 1 = “Yes”, 2 = “No”, etc.)

Will now look into the labels for the different levels that our factor variables can take:

attr(data$d70, "labels")

##       Very satisfied     Fairly satisfied   Not very satisfied 
##                    1                    2                    3 
## Not at all satisfied                   DK 
##                    4                    5

names(attr(data$d70, "labels"))

## [1] "Very satisfied"       "Fairly satisfied"     "Not very satisfied"  
## [4] "Not at all satisfied" "DK"

tibble(name_labels = names(attr(data$d70, "labels")),
       labels = attr(data$d70, "labels"))

## # A tibble: 5 × 2
##   name_labels          labels
##   <chr>                 <dbl>
## 1 Very satisfied            1
## 2 Fairly satisfied          2
## 3 Not very satisfied        3
## 4 Not at all satisfied      4
## 5 DK                        5

# Create a list of tibbles containing the labels and their associated name for each variable
list_label_tibbles <- 
  #Applies a function across all columns of a df and returns results as a list
  lapply(names(data), function(col_name) {
    labels <- attr(data[[col_name]], "labels")  # Extract labels
    name_labels <- names(labels)  # Extract label names
    # Create tibble with the extracted data only if labels exist
    if (!is.null(labels)) {
      tibble(name_labels = name_labels, labels = labels)} 
    else {NULL}  # Returns a NULL element for columns without labels
  })

# Giving to each element of the list as name the name of the variable
list_label_tibbles <- setNames(list_label_tibbles, names(data))

# For example
list_label_tibbles$d70

## # A tibble: 5 × 2
##   name_labels          labels
##   <chr>                 <dbl>
## 1 Very satisfied            1
## 2 Fairly satisfied          2
## 3 Not very satisfied        3
## 4 Not at all satisfied      4
## 5 DK                        5

Right column is what appears in our data (as a number). Left column is the label that we must assign to that number when we factorize

Initial cleaning and pre-processing

Recoding together Germany East and West because we are running analysis at the country level and they are coded separately.

unique(data$isocntry)

##  [1] "BE"   "DK"   "GR"   "ES"   "FI"   "FR"   "IE"   "IT"   "LU"   "NL"  
## [11] "AT"   "PT"   "SE"   "DE-W" "DE-E" "GB"   "BG"   "CY"   "CZ"   "EE"  
## [21] "HU"   "LV"   "LT"   "MT"   "PL"   "RO"   "SK"   "SI"   "HR"

data <- data |> 
  mutate(isocntry = case_when(
    isocntry %in% c("DE-W", "DE-E") ~ "DE",
    TRUE ~ isocntry))
# We might want to join the full name of the countries using the codelist df

Recode all factor NA’s

Separating variables for which we can use the attribute labels to factorize them and the variables for which this strategy cannot be used

data <- data |> 
  rename(friends_trans = sd1_7)

non_factor_variables <- c("serialid", "tnscntry", "isocntry", "d11", "q1_29", 
                          names(data)[startsWith(names(data), "sd1_")],
                          names(data)[startsWith(names(data), "sd2_")],
                          names(data)[startsWith(names(data), "qc5_")],
                          names(data)[startsWith(names(data), "qc12_")],
                          names(data)[startsWith(names(data), "qc13_")],
                          names(data)[startsWith(names(data), "qc15_")],
                          names(data)[startsWith(names(data), "d72_")],
                          "d8", "opls")
factor_variables <- setdiff(names(data), non_factor_variables)

Correctly encoding the factor variables that do not need any further cleaning

# Converting them to factors and assign them their labels automatically
data <- data |> 
  mutate(across(all_of(factor_variables), labelled::to_factor))

# Turning DK into NAs for all the factor variables
data <- data |> 
  mutate(across(all_of(factor_variables), ~ fct_na_level_to_value(., extra_levels = "DK")))

# Converting to numeric the variables that should be numeric
data <- data |> 
  mutate(age = as.numeric(d11),
         years_edu = as.numeric(d8)) |> 
  # Recoding correctly d8 education variable according to unique(data_raw$d11))
  mutate(years_edu = case_when(
    years_edu %in% c(0, 99) ~ NA, # Refusal and DK as NAs
    years_edu == 97 ~ 0, # No full time education = 0
    years_edu == 98 ~ age, # still studying = age
    TRUE ~ years_edu)) |> 
  select(-c(d11, d8))

# Converting q1_29 to factor without assigning labels (1 if non-Eu nationality, 0 if EU nationality)
# Same for sd2_
data <- data |> 
  mutate(nonEU_national = as.factor(q1_29),
         across(starts_with("sd2_"), ~ as.factor(.x))) |> 
  select(-q1_29)

Create new individual level measures

Dealing with the other variables on which we do some further pre-processing (feature engineering).

Creating a variable that counts the number of different minority groups a person has acquaintances with:

data <- data |>  
  mutate(across(starts_with("sd1"), ~ if_else(.x == 1, 1, 0))) |> 
  mutate(n_friends_minorities = sd1_1+sd1_2+sd1_3+sd1_4+sd1_5+sd1_6+sd1_8) |> 
  relocate(n_friends_minorities, .before="sd1_1") |> 
  select(-starts_with("sd1"))

Creating a variable that counts the number of actions against discrimination that you have taken in the last year:

data <- data |>  
  mutate(across(starts_with("qc5"), ~ if_else(.x == 1, 1, 0))) |> 
  mutate(n_actions_against_discri = qc5_1+qc5_2+qc5_3+qc5_4) |> 
  relocate(n_actions_against_discri, .after="qc3") |> 
  select(-starts_with("qc5"))

Building a discriminatory score:

data <- data |> 
  # Coding as NAs "it depends" and "DK"
  mutate(across(starts_with("qc12"), ~ if_else(.x >= 12, NA, .x))) |> 
  mutate(across(starts_with("qc13"), ~ if_else(.x >= 12, NA, .x))) |> 
  # Coding as 5 responses = indifferent
  mutate(across(starts_with("qc12"), ~ if_else(.x == 11, 5, .x))) |> 
  mutate(across(starts_with("qc13"), ~ if_else(.x == 11, 5, .x))) |> 
  # Modifying such that higher is more discriminatory
  mutate(roma_discri = 11 - rowMeans(cbind(qc12_1, qc13_1), na.rm = TRUE),
         black_discri = 11 - rowMeans(cbind(qc12_2, qc13_2), na.rm = TRUE),
         asian_discri = 11 - rowMeans(cbind(qc12_3, qc13_3), na.rm = TRUE),
         white_discri = 11 - rowMeans(cbind(qc12_4, qc13_4), na.rm = TRUE),
         jewish_discri = 11 - rowMeans(cbind(qc12_5, qc13_5), na.rm = TRUE),
         muslim_discri = 11 - rowMeans(cbind(qc12_6, qc13_6), na.rm = TRUE),
         buddihst_discri = 11 - rowMeans(cbind(qc12_7, qc13_7), na.rm = TRUE),
         christian_discri = 11 - rowMeans(cbind(qc12_8, qc13_8), na.rm = TRUE),
         atheist_discri = 11 - rowMeans(cbind(qc12_9, qc13_9), na.rm = TRUE),
         lgb_discri = 11 - rowMeans(cbind(qc12_10, qc13_10), na.rm = TRUE),
         trans_discri = 11 - rowMeans(cbind(qc12_11, qc13_11), na.rm = TRUE),
         intersex_discri = 11 - rowMeans(cbind(qc12_12, qc13_12), na.rm = TRUE),
         disability_discri = 11 - rowMeans(cbind(qc12_13, qc13_13), na.rm = TRUE),
         young_discri = 11 - rowMeans(cbind(qc12_14, qc13_14), na.rm = TRUE),
         old_discri = 11 - rowMeans(cbind(qc12_15, qc13_15), na.rm = TRUE)) |> 
  select(-starts_with("qc12")) |> 
  select(-starts_with("qc13"))

We delete the discrimination index against young and old people. This is because one of the question is: how comfortable you would feel if one of your children was in a love relationship with a person from group x. We don’t believe these are relevant to our target or discrimination score.

data <- data |> 
  select(-c("old_discri", "young_discri"))

Building a score of how supportive or anti lbtq rights you are:

data <- data |> 
  mutate(across(starts_with("qc15"), ~ if_else(.x == 5, NA, .x))) |>
  mutate(antilgbtq_rights = round(rowMeans(cbind(qc15_1, qc15_2, qc15_3), na.rm = TRUE), 2)) |> 
  select(-starts_with("qc15"))
# Scale of 1 to 4, 1 = supportive, 4 = homophobic

Aggregating together questions on whether you think your voice counts:

data <- data |> 
  mutate(across(starts_with("d72"), ~ if_else(.x > 4, NA, .x))) |>
  mutate(social_alienation = rowMeans(cbind(d72_1, d72_2), na.rm = TRUE)) |> 
  select(-starts_with("d72"))
# The higher the more people think their voice does not matter

This should be our final selection of variables.

We still need to rename them appropriately and check that all the NAs are correctly encoded by looking at the summary.

Further cleaning

Rename columns appropriately

data <- data |> 
  rename(
    country = isocntry,
    life_sat = d70,
    ethnic_minority = sd2_1,
    skincolor_minority = sd2_2,
    religious_minority = sd2_3,
    roma_minority = sd2_4,
    sexual_minority = sd2_5,
    disability_minority = sd2_6,
    religion = sd3,
    disc = qc2_15,
    disc_where = qc3,
    disc_contact = qc10,
    trans_docs = qc19,
    gender_docs = qc20,
    left_right = d1r1,
    marital_status = d7r1,
    gender = d10,
    occupation = d15a_r2,
    community = d25,
    phone_access = d43t,
    bill_issues = d60,
    internet_use = netuse,
    social_class = d63
  ) 

Now, we will transform the variable disc, a variable storing whether you have been subject to discrimination, which was originally coded in a negative way (1 = “Not mentioned”, 2 = “No, you haven’t been discriminated against”), into a positive dummy variable to make its interpretation more straightforward. The variable is 1 if a person has been subject to discrimination:

data <- data |> 
  mutate(suffered_discr = as.factor(ifelse(disc == "Not mentioned", 1, 0))) |> 
  select(-disc) |> 
  relocate(suffered_discr, .before = disc_where)

Move the variables around to have a ordered dataframe, and assign labels to help interpretation:

# Relocating to have a ordered df
data <- data |> 
  relocate(c("age", "gender", "years_edu","community", "marital_status", "occupation", "social_class", "religion", "nonEU_national", "phone_access", "bill_issues", "internet_use"), .after = country) |> 
  relocate(c("left_right", "social_alienation"), .after = polintr) |> 
  relocate(c("friends_trans", "n_friends_minorities", "n_actions_against_discri"), .after = gender_docs)

# Assigning labels to columns which have a difficult meaning
attr(data$gender, "label") <- NULL
attr(data$nonEU_national, "label") <- "OWNS A NON-EU NATIONALITY"
attr(data$social_alienation, "label") <- "HIGHER -> THINK THEIR VOICE DOESN'T  MATTER"
attr(data$ethnic_minority, "label") <- "ARE YOU PART OF X MINORITY" 
attr(data$skincolor_minority, "label") <- "ARE YOU PART OF X MINORITY" 
attr(data$religious_minority, "label") <- "ARE YOU PART OF X MINORITY" 
attr(data$roma_minority, "label") <- "ARE YOU PART OF X MINORITY" 
attr(data$sexual_minority, "label") <- "ARE YOU PART OF X MINORITY" 
attr(data$disability_minority, "label") <- "ARE YOU PART OF X MINORITY" 
attr(data$disability_minority, "label") <- "ARE YOU PART OF X MINORITY" 
attr(data$disability_minority, "label") <- "ARE YOU PART OF X MINORITY"
attr(data$suffered_discr, "label") <- "HAVE YOU BEEN SUBJECT TO DISCRIMINATION"
attr(data$n_friends_minorities, "label") <- "YOU KNOW PEOPLE FROM # NUMBER OF DIFFERENT MINORITES"
attr(data$n_actions_against_discri, "label") <- "YOU HAVE DONE # NUMBER OF DIFFERENT ACTIONS TO FIGHT DISCRIMINATIONS"
attr(data$roma_discri, "label") <- "HOW DISCRIMINATORY ARE YOU AGAINST X"
attr(data$black_discri, "label") <- "HOW DISCRIMINATORY ARE YOU AGAINST X"
attr(data$asian_discri, "label") <- "HOW DISCRIMINATORY ARE YOU AGAINST X"
attr(data$white_discri, "label") <- "HOW DISCRIMINATORY ARE YOU AGAINST X"
attr(data$jewish_discri, "label") <- "HOW DISCRIMINATORY ARE YOU AGAINST X"
attr(data$muslim_discri, "label") <- "HOW DISCRIMINATORY ARE YOU AGAINST X"
attr(data$buddihst_discri, "label") <- "HOW DISCRIMINATORY ARE YOU AGAINST X"
attr(data$christian_discri, "label") <- "HOW DISCRIMINATORY ARE YOU AGAINST X"
attr(data$atheist_discri, "label") <- "HOW DISCRIMINATORY ARE YOU AGAINST X"
attr(data$atheist_discri, "label") <- "HOW DISCRIMINATORY ARE YOU AGAINST X"
attr(data$lgb_discri, "label") <- "HOW DISCRIMINATORY ARE YOU AGAINST X"
attr(data$trans_discri, "label") <- "HOW DISCRIMINATORY ARE YOU AGAINST X"
attr(data$intersex_discri, "label") <- "HOW DISCRIMINATORY ARE YOU AGAINST X"
attr(data$disability_discri, "label") <- "HOW DISCRIMINATORY ARE YOU AGAINST X"
attr(data$antilgbtq_rights, "label") <- "HIGHER -> THEY OPPOSE MORE RIGHTS TO LGBTQ"

This is our dataset so far:

summary(data)

##     serialid       country               age          gender     
##  Min.   :    1   Length:27438       Min.   :15.00   Man  :12492  
##  1st Qu.: 6860   Class :character   1st Qu.:37.00   Woman:14946  
##  Median :13720   Mode  :character   Median :53.00                
##  Mean   :13720                      Mean   :51.56                
##  3rd Qu.:20579                      3rd Qu.:66.00                
##  Max.   :27438                      Max.   :98.00                
##                                                                  
##    years_edu                          community    
##  Min.   : 0.00   Rural area or village     : 8776  
##  1st Qu.:17.00   Small or middle sized town:10767  
##  Median :18.00   Large town                : 7881  
##  Mean   :19.66   NA's                      :   14  
##  3rd Qu.:22.00                                     
##  Max.   :90.00                                     
##  NA's   :387                                       
##                                 marital_status 
##  (Re-)Married (1-4 in d7)              :14673  
##  Single living with partner (5-8 in d7): 3321  
##  Single (9-10 in d7)                   : 4314  
##  Divorced or separated (11-12 in d7)   : 2243  
##  Widow (13-14 in d7)                   : 2712  
##  Other (SPONT.)                        :  107  
##  Refusal (SPONT.)                      :   68  
##                                   occupation  
##  Retired (4 in d15a)                   :8791  
##  Manual workers (15 to 18 in d15a)     :5883  
##  Other white collars (13 or 14 in d15a):3536  
##  Managers (10 to 12 in d15a)           :2911  
##  Self-employed (5 to 9 in d15a)        :1979  
##  Students (2 in d15a)                  :1676  
##  (Other)                               :2662  
##                             social_class                       religion    
##  The middle class of society      :13068   Catholic                :11198  
##  The working class of society     : 7233   Orthodox Christian      : 4016  
##  The lower middle class of society: 4070   Non believer or agnostic: 3694  
##  The upper middle class of society: 1889   Protestant              : 3031  
##  None (SPONTANEOUS)               :  229   Atheist                 : 2109  
##  (Other)                          :  389   (Other)                 : 3220  
##  NA's                             :  560   NA's                    :  170  
##  nonEU_national            phone_access               bill_issues   
##  0:27316        Mobile only      :14555   Most of the time  : 2054  
##  1:  122        Landline only    :  802   From time to time : 6538  
##                 Landline & mobile:11576   Almost never/never:18467  
##                 No telephone     :  505   Refusal (SPONT.)  :  379  
##                                                                     
##                                                                     
##                                                                     
##                      internet_use                   life_sat    
##  Everyday/Almost everyday  :19900   Very satisfied      : 7242  
##  Two or three times a week : 1701   Fairly satisfied    :15356  
##  About once a week         :  434   Not very satisfied  : 3736  
##  Two or three times a month:  164   Not at all satisfied: 1002  
##  Less often                :  350   NA's                :  102  
##  Never/No access           : 4311                               
##  No Internet access at all :  578                               
##        polintr               left_right   social_alienation ethnic_minority
##  Strong    : 4644   (1 - 4) Left  :7082   Min.   :1.000     0:26603        
##  Medium    :13701   (5 - 6) Centre:9313   1st Qu.:1.500     1:  835        
##  Low       : 4387   (7 -10) Right :6354   Median :2.000                    
##  Not at all: 4706   DK/Refusal    :4689   Mean   :2.329                    
##                                           3rd Qu.:3.000                    
##                                           Max.   :4.000                    
##                                           NA's   :1000                     
##  skincolor_minority religious_minority roma_minority sexual_minority
##  0:26923            0:26416            0:27001       0:27006        
##  1:  515            1: 1022            1:  437       1:  432        
##                                                                     
##                                                                     
##                                                                     
##                                                                     
##                                                                     
##  disability_minority suffered_discr
##  0:26756             0:23150       
##  1:  682             1: 4288       
##                                    
##                                    
##                                    
##                                    
##                                    
##                                                           disc_where   
##  In a public space                                             :  893  
##  At work                                                       :  788  
##  When looking for a job                                        :  642  
##  At a café, restaurant, bar or nightclub                       :  322  
##  By healthcare personnel (e.g. a receptionist, nurse or doctor):  308  
##  (Other)                                                       : 1054  
##  NA's                                                          :23431  
##                                                                       disc_contact 
##  The police                                                                 :8769  
##  A friend or family member                                                  :5386  
##  An equalities body or ombudsman (SPECIFY THE NAME ACCORDING TO THE COUNTRY):3795  
##  A lawyer                                                                   :2329  
##  Courts                                                                     :1225  
##  (Other)                                                                    :3601  
##  NA's                                                                       :2333  
##  trans_docs   gender_docs                friends_trans   n_friends_minorities
##  Yes :14463   Yes :10856   Yes                  : 2644   Min.   :0.000       
##  No  : 9695   No  :13395   No                   :23520   1st Qu.:1.000       
##  NA's: 3280   NA's: 3187   Refusal (SPONTANEOUS):  306   Median :3.000       
##                            NA's                 :  968   Mean   :3.029       
##                                                          3rd Qu.:5.000       
##                                                          Max.   :7.000       
##                                                                              
##  n_actions_against_discri  roma_discri      black_discri     asian_discri   
##  Min.   :0.0000           Min.   : 1.000   Min.   : 1.000   Min.   : 1.000  
##  1st Qu.:0.0000           1st Qu.: 1.500   1st Qu.: 1.000   1st Qu.: 1.000  
##  Median :0.0000           Median : 4.500   Median : 3.000   Median : 2.500  
##  Mean   :0.3855           Mean   : 4.631   Mean   : 3.579   Mean   : 3.387  
##  3rd Qu.:0.0000           3rd Qu.: 7.000   3rd Qu.: 5.500   3rd Qu.: 5.500  
##  Max.   :4.0000           Max.   :10.000   Max.   :10.000   Max.   :10.000  
##                           NA's   :757      NA's   :542      NA's   :592     
##   white_discri    jewish_discri   muslim_discri    buddihst_discri 
##  Min.   : 1.000   Min.   : 1.00   Min.   : 1.000   Min.   : 1.000  
##  1st Qu.: 1.000   1st Qu.: 1.00   1st Qu.: 1.000   1st Qu.: 1.000  
##  Median : 1.000   Median : 2.00   Median : 4.000   Median : 3.000  
##  Mean   : 1.785   Mean   : 3.18   Mean   : 4.363   Mean   : 3.576  
##  3rd Qu.: 2.000   3rd Qu.: 5.00   3rd Qu.: 6.500   3rd Qu.: 5.500  
##  Max.   :10.000   Max.   :10.00   Max.   :10.000   Max.   :10.000  
##  NA's   :403      NA's   :603     NA's   :662      NA's   :761     
##  christian_discri atheist_discri     lgb_discri      trans_discri   
##  Min.   : 1.000   Min.   : 1.000   Min.   : 1.000   Min.   : 1.000  
##  1st Qu.: 1.000   1st Qu.: 1.000   1st Qu.: 1.000   1st Qu.: 2.000  
##  Median : 1.000   Median : 1.500   Median : 4.000   Median : 5.000  
##  Mean   : 1.904   Mean   : 2.802   Mean   : 4.399   Mean   : 4.934  
##  3rd Qu.: 2.000   3rd Qu.: 4.000   3rd Qu.: 6.500   3rd Qu.: 7.500  
##  Max.   :10.000   Max.   :10.000   Max.   :10.000   Max.   :10.000  
##  NA's   :414      NA's   :517      NA's   :536      NA's   :1081    
##  intersex_discri  disability_discri antilgbtq_rights
##  Min.   : 1.000   Min.   : 1.000    Min.   :1.000   
##  1st Qu.: 1.500   1st Qu.: 1.000    1st Qu.:1.000   
##  Median : 5.000   Median : 2.500    Median :2.000   
##  Mean   : 4.825   Mean   : 2.977    Mean   :2.147   
##  3rd Qu.: 7.500   3rd Qu.: 4.500    3rd Qu.:3.000   
##  Max.   :10.000   Max.   :10.000    Max.   :4.000   
##  NA's   :1331     NA's   :585       NA's   :791

Our dataset seems quite balanced across the different countries

data |> count(country)

## # A tibble: 28 × 2
##    country     n
##    <chr>   <int>
##  1 AT       1027
##  2 BE       1028
##  3 BG       1032
##  4 CY        503
##  5 CZ       1008
##  6 DE       1537
##  7 DK       1004
##  8 EE       1003
##  9 ES       1005
## 10 FI       1003
## # ℹ 18 more rows

Although some levels of occupation have more observations than others, we have enough observations for each level so we do not need to aggregate over different levels for the variable occupation

data |> count(occupation)

## # A tibble: 8 × 2
##   occupation                                 n
##   <fct>                                  <int>
## 1 Self-employed (5 to 9 in d15a)          1979
## 2 Managers (10 to 12 in d15a)             2911
## 3 Other white collars (13 or 14 in d15a)  3536
## 4 Manual workers (15 to 18 in d15a)       5883
## 5 House persons (1 in d15a)               1358
## 6 Unemployed (3 in d15a)                  1304
## 7 Retired (4 in d15a)                     8791
## 8 Students (2 in d15a)                    1676

Very few people declare themselves to be the higher class of society, nonetheless we keep this level because it makes sense.

data |> count(social_class)

## # A tibble: 9 × 2
##   social_class                          n
##   <fct>                             <int>
## 1 The working class of society       7233
## 2 The lower middle class of society  4070
## 3 The middle class of society       13068
## 4 The upper middle class of society  1889
## 5 The higher class of society         157
## 6 Other (SPONTANEOUS)                  59
## 7 None (SPONTANEOUS)                  229
## 8 Refusal (SPONTANEOUS)               173
## 9 <NA>                                560

Given the low number of observations in some categories of the variable religion and the high number of categories we opt to aggregate some of them, otherwise we risk too much noise in our models.

data |> count(religion)

## # A tibble: 16 × 2
##    religion                     n
##    <fct>                    <int>
##  1 Catholic                 11198
##  2 Orthodox Christian        4016
##  3 Protestant                3031
##  4 Other Christian           1183
##  5 Jewish                      58
##  6 Muslim - Shia               77
##  7 Muslim - Sunni             178
##  8 Other Muslim               137
##  9 Sikh                        13
## 10 Buddhist                    58
## 11 Hindu                       32
## 12 Atheist                   2109
## 13 Non believer or agnostic  3694
## 14 Other                     1171
## 15 Refusal (SPONTANEOUS)      313
## 16 <NA>                       170

# Checking if the means of groups we are going to aggregate are similar
data |> 
  mutate(trans_docs=as.numeric(trans_docs)) |> 
  summarise(mean = mean(trans_docs, na.rm=TRUE), .by = religion) |>
  arrange(mean)

## # A tibble: 16 × 2
##    religion                  mean
##    <fct>                    <dbl>
##  1 Sikh                      1.25
##  2 Atheist                   1.27
##  3 Non believer or agnostic  1.28
##  4 Protestant                1.30
##  5 Buddhist                  1.33
##  6 Other                     1.37
##  7 Jewish                    1.38
##  8 Other Christian           1.43
##  9 Catholic                  1.43
## 10 Hindu                     1.44
## 11 Muslim - Shia             1.45
## 12 <NA>                      1.45
## 13 Other Muslim              1.49
## 14 Refusal (SPONTANEOUS)     1.49
## 15 Muslim - Sunni            1.54
## 16 Orthodox Christian        1.58

We are going to group together atheist with agnostic. We are also going to put sikh, buddhists, jewish and hindu into the other category because there are only very few observations. Finally we are going to group together all muslims.

data <- data |> 
  mutate(religion = fct_collapse(religion,
                                "Non-believers" = c("Atheist", "Non believer or agnostic"),
                                "Other" = c("Sikh", "Buddhist", "Jewish", "Hindu", "Other"),
                                "Muslim" = c("Muslim - Shia", "Muslim - Sunni", "Other Muslim")))

#This is what we end up with
data |> count(religion)

## # A tibble: 9 × 2
##   religion                  n
##   <fct>                 <int>
## 1 Catholic              11198
## 2 Orthodox Christian     4016
## 3 Protestant             3031
## 4 Other Christian        1183
## 5 Other                  1332
## 6 Muslim                  392
## 7 Non-believers          5803
## 8 Refusal (SPONTANEOUS)   313
## 9 <NA>                    170

Recoding missing values

In this section we will convert all levels of factors that are redundant to NAs.

First we will substitute de “Refusal” answers in the remaining variables as NAs. If the respondent refuses to answer, we treat is as equivalent to not knowing his or her answer. We exclude friends_trans because for that category it could be worth analyzing refusals.

factor_variables <- names(data)[sapply(data, is.factor)]

data <- data %>%
  mutate(across(
    all_of(setdiff(factor_variables, "friends_trans")),  # Exclude "friends_trans"
    ~ {
      # Look for levels that contain "REFUSAL"
      refusal_levels <- grep("REFUSAL", levels(.), value = TRUE, ignore.case = TRUE)
      # If there are levels containing "REFUSAL"
      if (length(refusal_levels) > 0) { 
        # Convert in NA
        fct_recode(., NULL = refusal_levels)  
      } 
      # if not remain without changes
      else {
        .  
      }
    }
  ))

Mutate NONE level of social_class to NA.

We also turn marital_status’s level OTHER into NAs, as it is difficult to give it any other meaning. Same for social class

data <- data %>%
  mutate(social_class = fct_recode(social_class, NULL = "None (SPONTANEOUS)"))

data <- data %>%
  mutate(marital_status = fct_recode(marital_status, NULL = "Other (SPONT.)"),
         social_class = fct_recode(social_class, NULL = "Other (SPONTANEOUS)"))

plot_intro(data)

plot_missing(data)

The total number of missing observations is quite low (4%) so missing data should not be a huge problem.

Given the high percentage of missing values we delete disc_where. This variable was expected to have a high percentage of missingness as it is a question that applies only to people that have been subject to discrimination. So it is not a variable of huge importance.

The variable with the second highest number of missing data is left_right, given the importance of this variable it is probably worth to impute those values.

Next there are the variables trans_docs and gender_docs which are respectively our target variable and a closely related variable

Then there is disc_contact which will be deleted as well, given that it does not add meaningful information to our analysis.

data <- data |> 
  select(-c("disc_where", "disc_contact"))

Sociodemographic variables

data_summary <- data |> 
  mutate(trans_docs = fct_na_value_to_level(trans_docs, "Don't Know")) |> 
  count(trans_docs, gender, name = "n") |> 
  group_by(gender) |>  
  mutate(percentage = n / sum(n))

ggplot(data_summary, aes(x = gender, y = percentage, fill = trans_docs)) +
  geom_bar(stat = "identity", position = "fill") +
  geom_text(aes(label = scales::percent(percentage, accuracy = 1)),  
            position = position_stack(vjust = 0.5), size = 4) + 
  scale_y_continuous(labels = percent_format()) + 
  scale_fill_manual(
    values = c("Yes" = "#2ECC71", "No" = "#E74C3C", "Don't Know" = "#3498DB")
  ) +
  labs(
    title = "Support for the right of trans people to change their gender in civil documents",
    x = "Gender",
    y = "Proportion of Responses",
    fill = "Support the right"
  ) +
  coord_flip() +  
  theme_minimal()+
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )

# Create contingency table
contingency_table <- table(data$gender, data$trans_docs)

# Calculate chisquare test
chisq_test <- chisq.test(contingency_table)
print(chisq_test)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  contingency_table
## X-squared = 92.715, df = 1, p-value < 2.2e-16

data |> 
  mutate(trans_docs = fct_na_value_to_level(trans_docs, "Don't Know")) |>
  mutate(age_bin = cut(age, breaks = seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE), by = 10), include.lowest = TRUE)) |> 
  # There are 4 observations above 95 y.o. that do not fall within any of the previously difined bins. Given their low number we drop them
  drop_na(age_bin) |> 
  count(age_bin, trans_docs) |> 
  group_by(age_bin) |> 
  mutate(percentage = n / sum(n) * 100) |>  
  ggplot(aes(x = age_bin, y = percentage, fill = trans_docs)) +
  geom_bar(stat = "identity", position = "stack") +
  coord_flip() +
  labs(x = "Age Bins", y = "Proportion of responses", fill = "Support the right",
       title = "Support for the right of trans people to change their gender in civil documents") +
  scale_fill_manual(
    values = c("Yes" = "#2ECC71", "No" = "#E74C3C", "Don't Know" = "#3498DB")
  ) +
  theme_minimal() +
   theme(
    panel.grid = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )

data |> 
  drop_na(trans_docs) |> 
  mutate(age_bin = cut(age, breaks = seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE), by = 10), include.lowest = TRUE)) |> 
  # There are 4 observations above 95 y.o. that do not fall within any of the previously defined bins. Given their low number we drop them
  drop_na(age_bin) |> 
  count(age_bin, trans_docs) |> 
  group_by(age_bin) |> 
  mutate(percentage = n / sum(n) * 100) |>  
  ggplot(aes(x = age_bin, y = percentage, fill = trans_docs)) +
  geom_bar(stat = "identity", position = "stack") +
  coord_flip() +
  labs(x = "Age Bins", y = "Proportion of responses", fill = "Support the right",
       title = "Support for the right of trans people to change their gender in civil documents") +
  scale_fill_manual(
    values = c("Yes" = "#2ECC71", "No" = "#E74C3C", "Don't Know" = "#3498DB")
  ) +
  theme_minimal() +
   theme(
    panel.grid = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )

data_summary <- data |> 
  mutate(trans_docs = fct_na_value_to_level(trans_docs, "Don't Know")) |>
  mutate(religious_minority = fct_recode(religious_minority, "Yes" = "1", "No" = "0")) |> 
  count(trans_docs, religion, name = "n") |> 
  group_by(religion) |>  
  mutate(percentage = n / sum(n))

ggplot(data_summary, aes(x = reorder(religion, ifelse(trans_docs == "Yes", -percentage, 0)), y = percentage, fill = trans_docs)) +
  geom_bar(stat = "identity", position = "fill") +
  geom_text(aes(label = scales::percent(percentage, accuracy = 1)),  
            position = position_stack(vjust = 0.5), size = 4) + 
  scale_y_continuous(labels = percent_format()) + 
  scale_fill_manual(
    values = c("Yes" = "#2ECC71", "No" = "#E74C3C", "Don't Know" = "#3498DB")
  ) +
  labs(
    title = "Support for the right of trans people to change their gender in civil documents",
    x = "Are you part of a religious minority?",
    y = "Proportion of Responses",
    fill = "Support the right"
  ) +
  coord_flip() +  
  theme_minimal()+
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )

# Create contingency table
contingency_table <- table(data$religion, data$trans_docs)

# Calculate chisquare test
chisq_test <- chisq.test(contingency_table)
print(chisq_test)

## 
##  Pearson's Chi-squared test
## 
## data:  contingency_table
## X-squared = 973.3, df = 6, p-value < 2.2e-16

data_summary <- data |> 
  mutate(trans_docs = fct_na_value_to_level(trans_docs, "Don't Know")) |>
  mutate(religious_minority = fct_recode(religious_minority, "Yes" = "1", "No" = "0")) |> 
  count(trans_docs, religious_minority, name = "n") |> 
  group_by(religious_minority) |>  
  mutate(percentage = n / sum(n))  

ggplot(data_summary, aes(x = religious_minority, y = percentage, fill = trans_docs)) +
  geom_bar(stat = "identity", position = "fill") +
  geom_text(aes(label = scales::percent(percentage, accuracy = 1)),  
            position = position_stack(vjust = 0.5), size = 4) + 
  scale_y_continuous(labels = percent_format()) + 
  scale_fill_manual(
    values = c("Yes" = "#2ECC71", "No" = "#E74C3C", "Don't Know" = "#3498DB")
  ) +
  labs(
    title = "Support for the right of trans people to change their gender in civil documents",
    x = "Are you part of a religious minority?",
    y = "Proportion of Responses",
    fill = "Support the right"
  ) +
  coord_flip() +  
  theme_minimal()+
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )

data_summary <- data |> 
  mutate(trans_docs = fct_na_value_to_level(trans_docs, "Don't Know")) |>
  count(trans_docs, left_right, name = "n") |> 
  group_by(left_right) |>  
  mutate(percentage = n / sum(n))  

ggplot(data_summary, aes(x = left_right, y = percentage, fill = trans_docs)) +
  geom_bar(stat = "identity", position = "fill") +
  geom_text(aes(label = scales::percent(percentage, accuracy = 1)),  
            position = position_stack(vjust = 0.5), size = 4) + 
  scale_y_continuous(labels = percent_format()) + 
  scale_fill_manual(
    values = c("Yes" = "#2ECC71", "No" = "#E74C3C", "Don't Know" = "#3498DB")
  ) +
  labs(
    title = "Support for the right of trans people to change their gender in civil documents",
    x = "Political affiliation",
    y = "Proportion of Responses",
    fill = "Support the right"
  ) +
  coord_flip() +  
  theme_minimal()+
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )

data_summary <- data |> 
  mutate(trans_docs = fct_na_value_to_level(trans_docs, "Don't Know")) |>
  count(trans_docs, social_class, name = "n") |> 
  group_by(social_class) |>  
  mutate(percentage = n / sum(n))  

ggplot(data_summary, aes(x = social_class, y = percentage, fill = trans_docs)) +
  geom_bar(stat = "identity", position = "fill") +
  geom_text(aes(label = scales::percent(percentage, accuracy = 1)),  
            position = position_stack(vjust = 0.5), size = 4) + 
  scale_y_continuous(labels = percent_format()) + 
  scale_fill_manual(
    values = c("Yes" = "#2ECC71", "No" = "#E74C3C", "Don't Know" = "#3498DB")
  ) +
  labs(
    title = "Support for the right of trans people to change \ntheir gender in civil documents",
    x = "You identify as:",
    y = "Proportion of Responses",
    fill = "Support the right"
  ) +
  coord_flip() +  
  theme_minimal()+
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )

# limit to main age ranges in age completed
data |> 
  mutate(main_educ = ifelse(between(years_edu, 15, 30), 1, 0)) |> 
  group_by(main_educ) |> 
  summarise(count = n(),
            prop = count/nrow(data)*100)

## # A tibble: 3 × 3
##   main_educ count  prop
##       <dbl> <int> <dbl>
## 1         0  3188 11.6 
## 2         1 23863 87.0 
## 3        NA   387  1.41

data |> 
  filter(between(years_edu, 15, 30)) |> 
  mutate(trans_docs = fct_na_value_to_level(trans_docs, "Don't Know")) |>
  count(trans_docs, years_edu, name = "n") |> 
  group_by(years_edu) |>  
  mutate(percentage = n / sum(n)) |> 
  ggplot(aes(x = years_edu, y = percentage, fill = trans_docs)) +
      geom_bar(stat = "identity", position = "fill") +
    geom_text(aes(label = scales::percent(percentage, accuracy = 1)),  
            position = position_stack(vjust = 0.5), size = 3) + 
  scale_y_continuous(labels = percent_format()) + 
  scale_fill_manual(
    values = c("Yes" = "#2ECC71", "No" = "#E74C3C", "Don't Know" = "#3498DB")
  ) +
  labs(
    title = "Distribution of support by education",
    x = "Age completed education, from 15 to 30 years old",
    y = "Proportion of Responses",
    fill = "Support the right"
  ) +
  theme_minimal()+
  theme(
    panel.grid = element_blank())

data_summary <- data |> 
  mutate(trans_docs = fct_na_value_to_level(trans_docs, "Don't Know")) |>
  count(trans_docs, community, name = "n") |> 
  group_by(community) |>  
  mutate(percentage = n / sum(n))  

ggplot(data_summary, aes(x = community, y = percentage, fill = trans_docs)) +
  geom_bar(stat = "identity", position = "fill") +
  geom_text(aes(label = scales::percent(percentage, accuracy = 1)),  
            position = position_stack(vjust = 0.5), size = 4) + 
  scale_y_continuous(labels = percent_format()) + 
  scale_fill_manual(
    values = c("Yes" = "#2ECC71", "No" = "#E74C3C", "Don't Know" = "#3498DB")
  ) +
  labs(
    title = "Support for the right of trans people to change their \ngender in civil documents",
    x = "Place of residence",
    y = "Proportion of Responses",
    fill = "Support the right"
  ) +
  coord_flip() +  
  theme_minimal()+
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )

Other individual-level variables

data_summary <- data |> 
  mutate(trans_docs = fct_na_value_to_level(trans_docs, "Don't Know")) |>
  count(trans_docs, friends_trans, name = "n") |> 
  group_by(friends_trans) |>  
  mutate(percentage = n / sum(n))  

ggplot(data_summary, aes(x = friends_trans, y = percentage, fill = trans_docs)) +
  geom_bar(stat = "identity", position = "fill") +
  geom_text(aes(label = scales::percent(percentage, accuracy = 1)),  
            position = position_stack(vjust = 0.5), size = 4) + 
  scale_y_continuous(labels = percent_format()) + 
  scale_fill_manual(
    values = c("Yes" = "#2ECC71", "No" = "#E74C3C", "Don't Know" = "#3498DB")
  ) +
  labs(
    title = "Support for the right of trans people to change their \ngender in civil documents",
    x = "Do you know a transgender person?",
    y = "Proportion of Responses",
    fill = "Support the right"
  ) +
  coord_flip() +  
  theme_minimal()+
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )

Number of persons present during interview

library(gmodels)

# Convert NA in a explicit category
merged_paradata$trans_docs <- as.character(merged_paradata$trans_docs)   
merged_paradata$trans_docs[is.na(merged_paradata$trans_docs)] <- "DK"  
merged_paradata$trans_docs <- as.factor(merged_paradata$trans_docs)  

# CrossTable
CrossTable(merged_paradata$p4, merged_paradata$trans_docs,
           digits=2, 
           expected=F, 
           asresid=T, 
           chisq=TRUE, 
           prop.chisq=F, 
           format="SPSS")

## 
##    Cell Contents
## |-------------------------|
## |                   Count |
## |             Row Percent |
## |          Column Percent |
## |           Total Percent |
## |           Adj Std Resid |
## |-------------------------|
## 
## Total Observations in Table:  27438 
## 
##                                  | merged_paradata$trans_docs 
##               merged_paradata$p4 |       DK  |       No  |      Yes  | Row Total | 
## ---------------------------------|-----------|-----------|-----------|-----------|
## Two (interviewer and respondent) |     2806  |     8146  |    12428  |    23380  | 
##                                  |    12.00% |    34.84% |    53.16% |    85.21% | 
##                                  |    85.55% |    84.02% |    85.93% |           | 
##                                  |    10.23% |    29.69% |    45.29% |           | 
##                                  |     0.58  |    -4.10  |     3.54  |           | 
## ---------------------------------|-----------|-----------|-----------|-----------|
##                            Three |      428  |     1316  |     1758  |     3502  | 
##                                  |    12.22% |    37.58% |    50.20% |    12.76% | 
##                                  |    13.05% |    13.57% |    12.16% |           | 
##                                  |     1.56% |     4.80% |     6.41% |           | 
##                                  |     0.52  |     2.97  |    -3.19  |           | 
## ---------------------------------|-----------|-----------|-----------|-----------|
##                             Four |       36  |      189  |      218  |      443  | 
##                                  |     8.13% |    42.66% |    49.21% |     1.61% | 
##                                  |     1.10% |     1.95% |     1.51% |           | 
##                                  |     0.13% |     0.69% |     0.79% |           | 
##                                  |    -2.50  |     3.25  |    -1.49  |           | 
## ---------------------------------|-----------|-----------|-----------|-----------|
##                     Five or more |       10  |       44  |       59  |      113  | 
##                                  |     8.85% |    38.94% |    52.21% |     0.41% | 
##                                  |     0.30% |     0.45% |     0.41% |           | 
##                                  |     0.04% |     0.16% |     0.22% |           | 
##                                  |    -1.02  |     0.80  |    -0.11  |           | 
## ---------------------------------|-----------|-----------|-----------|-----------|
##                     Column Total |     3280  |     9695  |    14463  |    27438  | 
##                                  |    11.95% |    35.33% |    52.71% |           | 
## ---------------------------------|-----------|-----------|-----------|-----------|
## 
##  
## Statistics for All Table Factors
## 
## 
## Pearson's Chi-squared test 
## ------------------------------------------------------------
## Chi^2 =  26.44717     d.f. =  6     p =  0.0001837381 
## 
## 
##  
##        Minimum expected frequency: 13.50827

The table shows a significant association between the number of people present during the interview and support for transgender people to change their documents. The actual support does not vary too much though, they are all within 3% of each other.

Duration of interview

Regarding the duration of the interview we have 2 options p3 and p3r. We will use p3r to have a simpler analysis

# CrossTable
CrossTable(merged_paradata$p3r, merged_paradata$trans_docs,
           digits=2, 
           expected=F, 
           asresid=T, 
           chisq=TRUE, 
           prop.chisq=F, 
           format="SPSS")

## 
##    Cell Contents
## |-------------------------|
## |                   Count |
## |             Row Percent |
## |          Column Percent |
## |           Total Percent |
## |           Adj Std Resid |
## |-------------------------|
## 
## Total Observations in Table:  27438 
## 
##                     | merged_paradata$trans_docs 
## merged_paradata$p3r |       DK  |       No  |      Yes  | Row Total | 
## --------------------|-----------|-----------|-----------|-----------|
##    Up to 14 minutes |       73  |      224  |      260  |      557  | 
##                     |    13.11% |    40.22% |    46.68% |     2.03% | 
##                     |     2.23% |     2.31% |     1.80% |           | 
##                     |     0.27% |     0.82% |     0.95% |           | 
##                     |     0.85  |     2.43  |    -2.88  |           | 
## --------------------|-----------|-----------|-----------|-----------|
##     15 - 29 minutes |      782  |     2284  |     3456  |     6522  | 
##                     |    11.99% |    35.02% |    52.99% |    23.77% | 
##                     |    23.84% |    23.56% |    23.90% |           | 
##                     |     2.85% |     8.32% |    12.60% |           | 
##                     |     0.10  |    -0.61  |     0.52  |           | 
## --------------------|-----------|-----------|-----------|-----------|
##     30 - 44 minutes |     1429  |     4388  |     6305  |    12122  | 
##                     |    11.79% |    36.20% |    52.01% |    44.18% | 
##                     |    43.57% |    45.26% |    43.59% |           | 
##                     |     5.21% |    15.99% |    22.98% |           | 
##                     |    -0.75  |     2.66  |    -2.06  |           | 
## --------------------|-----------|-----------|-----------|-----------|
##     45 - 59 minutes |      703  |     1993  |     3039  |     5735  | 
##                     |    12.26% |    34.75% |    52.99% |    20.90% | 
##                     |    21.43% |    20.56% |    21.01% |           | 
##                     |     2.56% |     7.26% |    11.08% |           | 
##                     |     0.80  |    -1.04  |     0.48  |           | 
## --------------------|-----------|-----------|-----------|-----------|
##     60 - 74 minutes |      202  |      529  |      908  |     1639  | 
##                     |    12.32% |    32.28% |    55.40% |     5.97% | 
##                     |     6.16% |     5.46% |     6.28% |           | 
##                     |     0.74% |     1.93% |     3.31% |           | 
##                     |     0.48  |    -2.67  |     2.25  |           | 
## --------------------|-----------|-----------|-----------|-----------|
##     75 - 89 minutes |       49  |      123  |      276  |      448  | 
##                     |    10.94% |    27.46% |    61.61% |     1.63% | 
##                     |     1.49% |     1.27% |     1.91% |           | 
##                     |     0.18% |     0.45% |     1.01% |           | 
##                     |    -0.67  |    -3.52  |     3.80  |           | 
## --------------------|-----------|-----------|-----------|-----------|
## 90 minutes and more |       42  |      154  |      219  |      415  | 
##                     |    10.12% |    37.11% |    52.77% |     1.51% | 
##                     |     1.28% |     1.59% |     1.51% |           | 
##                     |     0.15% |     0.56% |     0.80% |           | 
##                     |    -1.16  |     0.76  |     0.02  |           | 
## --------------------|-----------|-----------|-----------|-----------|
##        Column Total |     3280  |     9695  |    14463  |    27438  | 
##                     |    11.95% |    35.33% |    52.71% |           | 
## --------------------|-----------|-----------|-----------|-----------|
## 
##  
## Statistics for All Table Factors
## 
## 
## Pearson's Chi-squared test 
## ------------------------------------------------------------
## Chi^2 =  36.84302     d.f. =  12     p =  0.0002368813 
## 
## 
##  
##        Minimum expected frequency: 49.61003

It appears that the longer the interview, the more likely to support the people are. The highest support being 75-89 minutes while the highest opposition and and highest DK were in <15 minute interviews.

Respondent cooperation

data_summary <- merged_paradata |> 
  mutate(trans_docs = fct_recode(trans_docs, "Don't Know" = "DK"))|>
  count(p5, trans_docs, name = "n") |> 
  group_by(p5) |>  
  mutate(percentage = n / sum(n))  

ggplot(data_summary, aes(x = p5, y = percentage, fill = trans_docs)) +
  geom_bar(stat = "identity", position = "fill") +
  geom_text(aes(label = scales::percent(percentage, accuracy = 1)),  
            position = position_stack(vjust = 0.5), size = 4) + 
  scale_y_continuous(labels = percent_format()) + 
  scale_fill_manual(
    values = c("Yes" = "#2ECC71", "No" = "#E74C3C", "Don't Know" = "#3498DB")
  ) +
  labs(
    title = "Support for the right of trans people to change their gender in civil documents",
    x = "Respondent cooperation in the interview",
    y = "Proportion of Responses",
    fill = "Support the right"
  ) +
  coord_flip() +  
  theme_minimal()+
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )

Respondents who engaged better in the survey were more likely to express support. This is the most clear paradata relationship.

Low cooperation had higher DK response and opposition to trans documents support.

Check correlation with DK of target

Run model to test what is correlated with DK response, use Cramers V for association between the factor variables and correlation for the numeric variables in the cleaned data.

1. Model the data for DK responses against the factor variables.

factor_subset <- dk_target %>%
  dplyr::select(target_NA, where(is.factor)) |> 
  mutate(target_NA = factor(target_NA))

# Run cramers V for factors

target_var <- "target_NA"
factor_variables <- names(factor_subset)
factor_variables <- factor_variables[factor_variables != target_var]

cramers_v_results <- list()

for (variable in factor_variables) {
    contingency_table <- table(factor_subset[[target_var]], 
                               factor_subset[[variable]])
    cramers_v <- CramerV(contingency_table)
    cramers_v_results[[variable]] <- cramers_v
}

# Create a Tibble for Results
cat_results <- tibble(variable = names(cramers_v_results),
       cramers_v = unlist(cramers_v_results)) |> 
  arrange(desc(cramers_v))
cat_results

## # A tibble: 22 × 2
##    variable       cramers_v
##    <chr>              <dbl>
##  1 internet_use      0.136 
##  2 social_class      0.0975
##  3 occupation        0.0852
##  4 marital_status    0.0831
##  5 gender_docs       0.0820
##  6 religion          0.0763
##  7 phone_access      0.0711
##  8 life_sat          0.0642
##  9 polintr           0.0621
## 10 friends_trans     0.0581
## # ℹ 12 more rows

This shows all of the factors have quite low association with the NA responses. The highest being internet use. This is a good sign that the missing values are random.

2. Look at the numeric variables and test correlations:

numeric_subset <- dk_target %>%
  dplyr::select(target_NA, where(is.numeric)) 

target_var <- "target_NA"
numeric_variables <- names(numeric_subset)[names(numeric_subset) != target_var]

cor_results <- list()

for (variable in numeric_variables) {
    correlation <- cor(numeric_subset$target_NA, 
                       numeric_subset[[variable]],
                       use = "pairwise.complete.obs")
    cor_results[[variable]] <- correlation
}

# Create a Tibble for Results
cor_results <- tibble(variable = names(cor_results),
       cor = unlist(cor_results)) |> 
  arrange(cor)
cor_results

## # A tibble: 7 × 2
##   variable                     cor
##   <chr>                      <dbl>
## 1 n_actions_against_discri -0.0988
## 2 n_friends_minorities     -0.0965
## 3 years_edu                -0.0544
## 4 social_alienation         0.0566
## 5 minority_discri           0.0635
## 6 age                       0.0982
## 7 antilgbtq_rights          0.108

Again, this is not too much cause for concern. We only have around 10% correlations, positive and negative with the target variable.

Modelling relationship with key variables and non-response

Here we run a logistic regression to test for the most important/significant variables. First, we make a subset of data with only our variables that were most correlated in the seciton above (use above +/- 9.5% correlation and above 10% cramers V) and also include gender and country.

cor_results |> filter(cor > 0.09 | cor < -0.9) |> pull(variable)

## [1] "age"              "antilgbtq_rights"

cat_results |> filter(cramers_v >0.1) |> pull(variable)

## [1] "internet_use"

Now model with those variables to see if they are significant in explaining the DK values. Gender will also be included as a key variable.

We run a stepwise AIC on these most correlated variables and then use the best model to test the overall model fit.

# first, run stepwise to determine most important 
# run base model 
dk_fit_null <- glm(target_NA ~ 1,
              data = dk_target, 
              family = "binomial")
# run full fit model of most correlated vars with country (also age^2 for good practice)
dk_fit <- glm(target_NA ~ gender + age + I(age*age) + antilgbtq_rights + internet_use + social_class + country_name, 
              data = dk_target, 
              family = "binomial")

# use stepwise to select best mode l
aic_1 <- MASS::stepAIC(dk_fit, scope = list(upper = dk_fit, 
                             lower = dk_fit_null),
        direction = "both", k = 2, trace=0) # forward based on AIC

# filter vars significant at 5% and calculate exponential
broom::tidy(aic_1) |> 
  filter(p.value<0.05) |> 
  mutate(oddsratio = exp(estimate)) |> 
  select(term, oddsratio, p.value)

## # A tibble: 29 × 3
##    term                                          oddsratio  p.value
##    <chr>                                             <dbl>    <dbl>
##  1 (Intercept)                                      0.0500 8.78e-40
##  2 genderWoman                                      1.20   2.01e- 5
##  3 I(age * age)                                     1.00   4.77e- 4
##  4 antilgbtq_rights                                 1.13   9.69e- 7
##  5 internet_useTwo or three times a week            1.21   2.72e- 2
##  6 internet_useTwo or three times a month           0.307  4.97e- 3
##  7 internet_useNever/No access                      1.28   2.49e- 4
##  8 internet_useNo Internet access at all            1.57   4.62e- 4
##  9 social_classThe lower middle class of society    0.766  4.62e- 5
## 10 social_classThe middle class of society          0.713  4.09e-11
## # ℹ 19 more rows

The best model is the one without age but includes age^2. So we run the full model to keep the lower level variables included.

The model has some interesting findings to who did not respond to our findings. Some things we can see with the odds ratios (excluding interpretation of countries as that is covered above):

• women are around 20% more likely to give NA responses

• Anti-lgbti rights people were about 10% more likely to not respond

• Internet use was significant, those who had not used the internet or did not have access, were less likely to respond compared to those who use the internet everyday.

• Social class was significant at all levels, this means the working class were more likely to NA respond than all other class levels.

Below is an option using undersampling (to reduce class imbalance of being a NA versus not being a NA for our target variable).

na_model <- glm(target_NA ~ .,
                      data = dk_target, 
                      family = "binomial")

broom::tidy(na_model) |> arrange(desc(estimate))

## # A tibble: 89 × 5
##    term                       estimate std.error statistic       p.value
##    <chr>                         <dbl>     <dbl>     <dbl>         <dbl>
##  1 country_nameBulgaria          1.26     0.346       3.63 0.000280     
##  2 country_nameSweden            0.822    0.320       2.56 0.0103       
##  3 religionMuslim                0.692    0.325       2.13 0.0335       
##  4 country_nameLatvia            0.687    0.321       2.14 0.0323       
##  5 country_nameCyprus            0.619    0.450       1.38 0.169        
##  6 country_nameGermany           0.552    0.303       1.82 0.0683       
##  7 country_nameUnited Kingdom    0.545    0.329       1.66 0.0977       
##  8 gender_docsNo                 0.536    0.0930      5.77 0.00000000796
##  9 country_nameFinland           0.507    0.328       1.54 0.123        
## 10 country_namePoland            0.501    0.326       1.54 0.124        
## # ℹ 79 more rows

Approach to test relationship to DK with splitting and undersampling:

# create new dataset with just vars we want, including age^2
dk_target2 <- dk_target |> 
  select(target_NA, age, antilgbtq_rights, internet_use, social_class, country_name) |> 
  mutate(age_sq = age*age, .after = age)

# to deal with some class imbalance, undersample from the majority group (NA = 0)
set.seed(123)
# Generate an index
index <- createDataPartition(dk_target2$target_NA, p = 0.7, list = FALSE, times = 1)

# Subset the dataframe
train <- dk_target2[index, ]
test <- dk_target2[-index, ]
# check splits 
prop.table(table(train$target_NA))

## 
##         0         1 
## 0.8787421 0.1212579

prop.table(table(test$target_NA))

## 
##         0         1 
## 0.8844612 0.1155388

Now undersample the majority

set.seed(123)
under <- ovun.sample(target_NA~., 
                     data=train, 
                     method = "under", 
                     N = 4000)$data 
    # limit the sample to 4000 obs as we have ~2300 for target minority class. 

table(under$target_NA)

## 
##    0    1 
## 2112 1888

#run a rf model to see how much of the NA values we can explain
rfunder <- randomForest(target_NA~., data=under)
rfunder

## 
## Call:
##  randomForest(formula = target_NA ~ ., data = under) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##           Mean of squared residuals: 0.239057
##                     % Var explained: 4.08

# check a simple logistic model oucomes for undersampled data
broom::tidy(glm(target_NA ~., 
                data = under, 
                family="binomial")) |> 
  filter(p.value<0.05) |> 
  mutate(oddsratio = exp(estimate)) |> 
  select(term, oddsratio, p.value)

## # A tibble: 15 × 3
##    term                                          oddsratio    p.value
##    <chr>                                             <dbl>      <dbl>
##  1 (Intercept)                                       0.451 0.0178    
##  2 antilgbtq_rights                                  1.13  0.00254   
##  3 internet_useTwo or three times a week             1.39  0.0173    
##  4 internet_useNever/No access                       1.55  0.000114  
##  5 internet_useNo Internet access at all             1.85  0.0120    
##  6 social_classThe middle class of society           0.681 0.00000455
##  7 social_classThe upper middle class of society     0.705 0.0297    
##  8 country_nameBelgium                               0.141 0.00000805
##  9 country_nameBulgaria                              2.62  0.0000802 
## 10 country_nameCyprus                                2.72  0.000779  
## 11 country_nameEstonia                               1.65  0.0454    
## 12 country_nameLatvia                                1.88  0.0125    
## 13 country_nameLithuania                             2.06  0.00589   
## 14 country_nameNetherlands                           0.423 0.00793   
## 15 country_namePoland                                2.23  0.00143

We see that only around 5% of the variance is explained through RF. However The undersampled data with a simplified dataframe finds that higher DK were associated with being anti lgbti rights, less frequent internet use and being working class. This is all useful to compare to overall response rates.

Testing the NA values against paradata

Here we extend the model to see if paradata is related to changed non-response.

para <- paradata |> 
  bind_cols(trans_docs = data$trans_docs,
            age = data$age,
            gender = data$gender,
            country = data$country) |> 
  mutate(target_NA = ifelse(is.na(trans_docs), 1, 0)) |> 
  select(-serialid, -trans_docs)

summary(para$target_NA)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.1195  0.0000  1.0000

str(para)

## tibble [27,438 × 9] (S3: tbl_df/tbl/data.frame)
##  $ p2       : Factor w/ 6 levels "Before 8 h","8 - 12 h",..: 4 4 3 4 2 4 5 2 4 4 ...
##   ..- attr(*, "label")= chr "TIME OF INTERVIEW"
##  $ p3       : Factor w/ 246 levels "2 minutes","3",..: 44 39 37 30 30 25 26 33 39 39 ...
##   ..- attr(*, "label")= chr "DURATION OF INTERVIEW"
##  $ p3r      : Factor w/ 7 levels "Up to 14 minutes",..: 4 3 3 3 3 2 2 3 3 3 ...
##   ..- attr(*, "label")= chr "DURATION OF INTERVIEW (RECODED)"
##  $ p4       : Factor w/ 5 levels "Two (interviewer and respondent)",..: 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "label")= chr "N OF PERSONS PRESENT DURING INTERVIEW"
##  $ p5       : Factor w/ 5 levels "Excellent","Fair",..: 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "label")= chr "RESPONDENT COOPERATION"
##  $ age      : num [1:27438] 51 62 38 29 63 41 48 88 44 45 ...
##  $ gender   : Factor w/ 2 levels "Man","Woman": 1 2 2 1 2 2 1 1 2 2 ...
##  $ country  : chr [1:27438] "BE" "BE" "BE" "BE" ...
##  $ target_NA: num [1:27438] 0 0 0 0 0 0 0 0 0 0 ...

First, again test correlation. We use Cramers V since we only have factors:

target_var <- "target_NA"
para_vars <- names(para)[names(para) != target_var]

cramers_v_results <- list()

for (variable in para_vars) {
    contingency_table <- table(para[[target_var]], 
                               para[[variable]])
    contingency_table_no_dk <- contingency_table[, colnames(contingency_table) != "DK"]
    cramers_v <- CramerV(contingency_table_no_dk) 
    # include this for no DK as it is all zero which creates errors
    cramers_v_results[[variable]] <- cramers_v
}


# Create a Tibble for Results
tibble(variable = names(cramers_v_results),
       cramers_v = unlist(cramers_v_results)) |> 
  arrange(desc(cramers_v))

## # A tibble: 8 × 2
##   variable cramers_v
##   <chr>        <dbl>
## 1 country     0.168 
## 2 age         0.123 
## 3 p5          0.106 
## 4 p3          0.0916
## 5 gender      0.0292
## 6 p4          0.0165
## 7 p2          0.0116
## 8 p3r         0.0113

Again, we see low correlation with the paradata and cramers V.

Second, we run another logistic model with just Paradata included.

na_para_model <- glm(target_NA ~ p2 + p3r + p4 + p5 + age + I(age*age) + gender + factor(country), 
                     data = para, 
                     family=binomial(link = "logit"))

broom::tidy(na_para_model) |> 
  filter(p.value<0.05) |> 
  mutate(oddsratio = exp(estimate)) |> 
  select(term, oddsratio, p.value)

## # A tibble: 25 × 3
##    term              oddsratio  p.value
##    <chr>                 <dbl>    <dbl>
##  1 (Intercept)          0.0439 4.65e-12
##  2 p4Four               0.694  4.20e- 2
##  3 p5Fair               1.42   2.05e-15
##  4 p5Average            1.74   8.91e-18
##  5 p5Bad                3.49   2.14e-22
##  6 age                  0.983  2.85e- 3
##  7 I(age * age)         1.00   3.43e- 9
##  8 genderWoman          1.14   1.07e- 3
##  9 factor(country)BE    0.179  6.77e- 9
## 10 factor(country)BG    5.30   2.16e-31
## # ℹ 15 more rows

We see that when controlling for age, gender and country, the survey cooperation variable is significant.

When compared to those with an excellent cooperation rating, people were more likely to respond. This increases with lacking cooperation. If cooperation was fair, there was around 40% more chance of NA, if cooperation was average, around 75% higher chance of NA. If cooperation was bad, there was about 350% higher chance of a NA response to our target variable. This variable also has risk of reverse causality, as higher NA respose is likely a reason to be rated as low cooperation.

Overall, this section found that people more likely to respond NA seems similar to those more likely to oppose trans rights to change papers. So there is some risk of under reporting of true support if people opposed are not responding. This is likely capturing some social desirability as people don’t want to appear discriminatory and hide their true views.

Interpeting results

To interpret the coefficients as odds ratios, anything above 1 indicates they are more likely to support the changes to civil documents for trans people.

The logistic regression shows that 26 statistically significant terms. This includes factors variables with multiple terms. Of our 29 scaled variables, the significant variables (at the 5% level) are:

More likely to OPPPOSE the right to change civil documents for trans people, include:

• men (compared to women)

• younger people (this is in a way unexpected)

• people who oppose legal rights to add a third-gender in official docs (compared to those who do)

• people who refused to answer whether they have trans friends (compared to those who do)

• people who are more financially sound (have fewer bill issues)

• people not satisfied with their life (compared to those who are ‘very satisfied’)

• people who identified as roma or an ethnic minority

• people who are more discriminatory against minorities

More likely to SUPPORT the the right

• older people

• women

• people who were self-employed (compared to all other occupation types incl students)

• unmarried (single) people

• non-believers (religious)

• people with a landline and mobile

• people who use the internet everyday/almost everyday (compared to all other internet use categories)

• people who reported being more left wing (compared to right wing)

• people who had friends in minority groups.

NUMERIC VARIABLES

Without scaling, we interpret logistic regression output as:

coef(simple_logistic): these coefficients represent the change in the log-odds for a one-unit increase in the corresponding independent variable

exp(coef(simple_logistic)): the exponential of the slope coefficient (exp(B)) tells us the change of the odds if the independent variable increases by one unit

If the data is scaled:

A coefficient B (from a logistic regression) now reflects the change in the log-odds for a one standard deviation increase in the predictor variable, rather than a one-unit increase in the original scale.

The odds ratio (exp(B)) now tells you how the odds change with a one standard deviation increase in the predictor variable, rather than a one-unit increase.

Here we compare the pre and post transformation of coefficients for interpretation:

# all the coefficients
head(coef(simple_logistic))

##                         (Intercept)                                 age 
##                          2.69303178                          0.20879798 
##                         genderWoman                           years_edu 
##                          0.26342657                         -0.01215087 
## communitySmall or middle sized town                 communityLarge town 
##                          0.01810837                         -0.04039684

# for their interpretation
head(exp(coef(simple_logistic)))

##                         (Intercept)                                 age 
##                          14.7764069                           1.2321960 
##                         genderWoman                           years_edu 
##                           1.3013817                           0.9879227 
## communitySmall or middle sized town                 communityLarge town 
##                           1.0182733                           0.9604082

FACTOR VARIABLES

Factor variables when exponentialised give the likelihood to respond no/yes relative to the base group in the factor variable. We use the function factors to view the levels and identify the base group.

To assist with interpretation of factor outputs:

Here is a function to identify all the base levels in each factor group. Useful to refer to during analysis of different logistic models.

extract_base_groups <- function(df) {

  # extract factor vars only and list to print values
  factor_vars <- names(final_data)[sapply(final_data, is.factor)]
  base_groups <- list()
  
  # loop through values to check all levels
  for (var in factor_vars) {
    if (length(levels(df[[var]])) > 0) { # Check step - should be TRUE for all selected. 
      contrast_matrix <- contrasts(df[[var]])
      if (is.null(contrast_matrix)) {
        base_groups[[var]] <- levels(df[[var]])[1] # check only - shouldn't ever run
      } else {
        base_group_level <- levels(df[[var]])[rowSums(abs(contrast_matrix)) == 0] # base group is the one with all zeroes, so we extract that
        base_groups[[var]] <- base_group_level
      }
    } else {
      base_groups[[var]] <- NA # Indicate non-factor or factor with no levels.
    }
  }
  return(base_groups)
}

# run function with 'final data' and print a tidy output
factor_bases <- extract_base_groups(final_data)

factor_bases <- as_tibble(factor_bases) |> 
  pivot_longer(cols = everything(), 
               names_to = "vars",
               values_to = "base_group")

print(factor_bases, n=Inf)

## # A tibble: 23 × 2
##    vars                base_group                    
##    <chr>               <chr>                         
##  1 gender              Man                           
##  2 community           Rural area or village         
##  3 marital_status      (Re-)Married (1-4 in d7)      
##  4 occupation          Self-employed (5 to 9 in d15a)
##  5 social_class        The working class of society  
##  6 religion            Catholic                      
##  7 nonEU_national      0                             
##  8 phone_access        Mobile only                   
##  9 bill_issues         Most of the time              
## 10 internet_use        Everyday/Almost everyday      
## 11 life_sat            Very satisfied                
## 12 polintr             Strong                        
## 13 left_right          (1 - 4) Left                  
## 14 ethnic_minority     0                             
## 15 skincolor_minority  0                             
## 16 religious_minority  0                             
## 17 roma_minority       0                             
## 18 sexual_minority     0                             
## 19 disability_minority 0                             
## 20 suffered_discr      0                             
## 21 gender_docs         Yes                           
## 22 friends_trans       Yes                           
## 23 trans_docs          Yes

Model performance - Simple logistic

Here, we test the models ability to classify people as supporting or not supporting the rights.

# AIC of model
simple_logistic$aic

## [1] 19544.18

#predicted probabilities of being a 1 (i.e. a refusal of possibility of trans doc)
predicted_probs <- simple_logistic$fitted.values
# same by running
# predicted_probs <- predict(simple_logistic, type = "response")
head(predicted_probs)

##         1         2         3         4         5         6 
## 0.3463520 0.4321645 0.5170989 0.5799662 0.9141149 0.5131253

# AUC of the model
auc(datalr$trans_docs, predicted_probs)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## Area under the curve: 0.8937

# Turning them to classes using custom threshold
predicted_classes <- ifelse(predicted_probs > 0.5, 1, 0)
head(predicted_classes)

## 1 2 3 4 5 6 
## 0 0 1 1 1 1

Get the confusion matrix:

cM <- confusionMatrix(as.factor(datalr$trans_docs), as.factor(predicted_classes))
cM

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction     0     1
##          0  7608  2087
##          1  2246 12217
##                                           
##                Accuracy : 0.8206          
##                  95% CI : (0.8157, 0.8255)
##     No Information Rate : 0.5921          
##     P-Value [Acc > NIR] : < 2e-16         
##                                           
##                   Kappa : 0.6277          
##                                           
##  Mcnemar's Test P-Value : 0.01638         
##                                           
##             Sensitivity : 0.7721          
##             Specificity : 0.8541          
##          Pos Pred Value : 0.7847          
##          Neg Pred Value : 0.8447          
##              Prevalence : 0.4079          
##          Detection Rate : 0.3149          
##    Detection Prevalence : 0.4013          
##       Balanced Accuracy : 0.8131          
##                                           
##        'Positive' Class : 0               
## 

The simple model predicted people’s response with 82% accuracy. That is quite decent for a basic first model. This shows the response is quite explainable with the data we have. The model is slightly better at predicting those who do not support than those who do. Seen by higher specificity than sensitivity (85% vs 77%). But they are both comparable which is good and confirms we don’t have great class imbalance.

Kappa is basically an accurancy score, but taking into account also the possibility of being right by chance. A value of 0.62 is considered quite good.

plot.roc(datalr$trans_docs, predicted_probs, 
         col="darkblue", 
         print.auc = TRUE,  
         auc.polygon=TRUE,
         max.auc.polygon=TRUE,
         auc.polygon.col="lightblue", 
         print.thres="best")

This suggests changing the threshold for assigning classes to 0.59

Anyway we are not really interested in the prediction power of this model, but we are more interested in understanding which variables are relevant or not

Final variable selection

From these hints above and general intuition, we will select the following variables as our individual level fixed effects in the mixed level model:

• age

• gender

• religion: given that the majority of people in our sample are catholic while the most important levels seem to be non_believers, muslims and Orthodox, instead of using all the levels we are going to use these 3 levels

• occupation: it is not clear which levels might be relevant so we will not include any of them

• Place you live does not seem to matter (hence no community variables in our final model)

• marital_status it is not clear which levels are important so we drop this variable

• social_class as our descriptive analysis was suggesting that this might be relevant and given that we are not using occupation, we will include all levels of social_class

• phone_access seems like it might be important, but its interpretation is difficult and I am not sure what this variable can be a proxy of, therefore we exclude it

• bill_issues will be recoded as a dummy using the level almost never. This means that those with 1s in our dummy experience difficult paying bills on a sporadic or constant basis. This can easily be interpreted as a proxy of whether a person is poor or not.

• internet_use on a similar concept as the previous point, this variable will be recoded as a dummy storing whether one uses internet every day or not.

• life_sat: people who are not satisfied in life (Not very satisfied & Not at all satisfied) seem to be really different from the rest. So a dummy will be created using those 2 levels

• left_right: right wing people seem significantly different than center or left wing so a dummy will be created for right wing people

• ethnic_minority and roma_minority seem quite important. Given that Roma is an ethnic minority we will keep only the former and discard the latter, which would be redundant information

• Whether you have suffered discrimination (suffered_discri) matters

• n_friends_minorites and friends_trans are very similar concepts. n_friends_minorites is derived summing the “yes” answers to “do you know a person that belong to this minority x”, but excluding the transgender minority whose answers are left as a standalone variable (friends_trans). We are going to use directly friends_trans as it is more directly related with our target, discarding n_friends_minorities that likely captures a similar effect

• antilgbtq_rights is clearly a good predictor. As gender_docs again captures a similar concept we are going to use only antilgbtq_rights

• minority_disci it is also closely related to the last point above, but we are going to keep it as a more general score of how discriminatory you are

And these will be our pre-defined fixed effects.

# Preparing the dataframe
complete_df <- complete_df |> 
  # Irrelevant variables
  select(-c(years_edu, community, marital_status, occupation,
            nonEU_national, phone_access, polintr, social_alienation,
            gender_docs, n_friends_minorities, 
            n_actions_against_discri, skincolor_minority,
            religious_minority, roma_minority, sexual_minority,
            disability_minority)) |> 
  # Religion
  mutate(non_believer = ifelse(religion == "Non-believers", 1, 0),
    muslim = ifelse(religion == "Muslim", 1, 0),
    orthodox = ifelse(religion == "Orthodox Christian", 1, 0)) |>
  select(-religion) |> 
  # Bill issues
  mutate(poor = ifelse(bill_issues != "Almost never/never", 1, 0)) |> 
  select(-bill_issues) |> 
  # Internet use
  mutate(everyday_internet = ifelse(internet_use == "Everyday/Almost everyday", 1, 0))  |> 
  select(-internet_use) |> 
  # Life satisfaction
  mutate(unhappy = ifelse(life_sat %in% c("Not very satisfied", "Not at all satisfied"), 1, 0)) |> 
  select(-life_sat) |> 
  # Ideology
  mutate(right_wing = ifelse(left_right == "(7 -10) Right", 1, 0)) |> 
  select(-left_right)

str(complete_df)

## 'data.frame':    27438 obs. of  21 variables:
##  $ country                : chr  "BE" "BE" "BE" "BE" ...
##  $ age                    : num  51 62 38 29 63 41 48 88 44 45 ...
##  $ gender                 : Factor w/ 2 levels "Man","Woman": 1 2 2 1 2 2 1 1 2 2 ...
##  $ social_class           : Factor w/ 5 levels "The working class of society",..: 2 3 2 4 4 3 3 3 3 4 ...
##  $ ethnic_minority        : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "label")= chr "ARE YOU PART OF X MINORITY"
##  $ suffered_discr         : Factor w/ 2 levels "0","1": 1 2 2 2 1 1 1 1 1 1 ...
##   ..- attr(*, "label")= chr "HAVE YOU BEEN SUBJECT TO DISCRIMINATION"
##  $ friends_trans          : Factor w/ 3 levels "Yes","No","Refusal (SPONTANEOUS)": 2 2 2 2 2 2 2 2 2 2 ...
##   ..- attr(*, "label")= chr "CONTACT: TRANSGENDER/TRANSSEXUAL"
##  $ antilgbtq_rights       : num  2 2.33 1.33 1.33 2.67 1.33 1.33 2.33 1.67 1 ...
##   ..- attr(*, "label")= chr "HIGHER -> THEY OPPOSE MORE RIGHTS TO LGBTQ"
##  $ minority_discri        : num  3.38 3.75 3.38 1 3.12 ...
##  $ trans_docs             : Factor w/ 2 levels "Yes","No": 1 1 1 1 1 1 1 2 1 1 ...
##   ..- attr(*, "label")= chr "TRANSGENDER - CHANGE CIVIL DOCUMENTS"
##  $ gdp_pc_ppp             : num  60452 60452 60452 60452 60452 ...
##  $ gender_inequality_index: num  0.048 0.048 0.048 0.048 0.048 0.048 0.048 0.048 0.048 0.048 ...
##  $ lgbt_policy_index      : num  10 10 10 10 10 10 10 10 10 10 ...
##  $ democracy_index        : num  76.4 76.4 76.4 76.4 76.4 76.4 76.4 76.4 76.4 76.4 ...
##  $ non_believer           : num  0 0 1 1 0 1 1 0 0 0 ...
##  $ muslim                 : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ orthodox               : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ poor                   : num  1 0 1 0 0 0 1 1 0 0 ...
##  $ everyday_internet      : num  1 1 1 1 1 1 1 0 1 1 ...
##  $ unhappy                : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ right_wing             : num  1 1 0 1 1 1 0 1 0 1 ...

Running the glmer models

To model, we will follow: Approach to multilevel model building based on Hox (2010)

1/ Null Model (Random Intercept only)

2/ Add independent Level 1 variables

3/ Add independent Level 2 variables

4/ Add random slopes

5/ (Cross-level) interactions

Each step, check whether your model is significantly improved compared to the previous one.

Adjusted ICC: The proportion of total variance explained by the grouping structure (random effects), adjusted to include the variance explained by the fixed effects

Unadjusted ICC: The proportion of total variance explained by the grouping structure (random effects), excluding any variance explained by the fixed effects in the model

complete_df <- complete_df |> 
  drop_na() |> # we have NAs only in our target
  mutate(trans_docs = as.numeric(trans_docs == "Yes"),
         ethnic_minority = as.numeric(ethnic_minority) - 1,
         suffered_discr = as.numeric(suffered_discr) - 1) 

# Scaling the variables is necessary otherwise the model cannot run
# However we loose interpretability 
scaled_df <- complete_df |>
  mutate(across(c(age, antilgbtq_rights, minority_discri, gdp_pc_ppp,
                  gender_inequality_index, lgbt_policy_index, 
                  democracy_index), ~scale(.x, center = TRUE, scale = TRUE)))

# Null model random intercept only
glmer_step1 <- glmer(trans_docs ~ (1|country),
                          data = scaled_df, 
                          family=binomial(link="logit"))
# Printing coefficients
nice_table(broom.mixed::tidy(glmer_step1))

# Random intercepts coefficients
ranef(glmer_step1)

## $country
##     (Intercept)
## AT  0.020782239
## BE  0.531337356
## BG -2.038432588
## CY -0.290724298
## CZ -0.658351324
## DE  0.875561535
## DK  0.910532764
## EE  0.009170019
## ES  1.763206743
## FI  0.443774385
## FR  0.870511332
## GB  0.402502267
## GR -0.051877595
## HR -0.631306153
## HU -2.093833734
## IE  0.554181536
## IT -0.539249991
## LT -0.741686697
## LU  0.854896968
## LV -0.510739696
## MT  1.282868516
## NL  1.406610447
## PL -0.447199501
## PT  0.782044709
## RO -1.599375278
## SE  0.500727881
## SI -0.334788326
## SK -1.293857996
## 
## with conditional variances for "country"

# AIC
AIC(glmer_step1)

## [1] 28212.65

# AUC
predicted_probs <- predict(glmer_step1)
auc(scaled_df$trans_docs, predicted_probs)

## Area under the curve: 0.7405

# ICC
performance::icc(glmer_step1)

## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.227
##   Unadjusted ICC: 0.227

# Random intercept + individual fixed effects
glmer_step2 <- glmer(trans_docs ~ age + gender + social_class +
                       ethnic_minority + suffered_discr + 
                       friends_trans + antilgbtq_rights +
                       minority_discri + non_believer + muslim +
                       orthodox + poor + everyday_internet + 
                       unhappy + right_wing + (1|country),
                     data = scaled_df,
                     family = binomial(link="logit"))
# Printing coefficients
nice_table(broom.mixed::tidy(glmer_step2))

# Random intercepts coefficients
ranef(glmer_step2)

## $country
##    (Intercept)
## AT -0.02120341
## BE  0.22958175
## BG -1.11582225
## CY  0.19511419
## CZ -0.48792517
## DE  0.37873507
## DK  0.05356456
## EE  0.58148184
## ES  1.10243054
## FI  0.13358429
## FR  0.29652703
## GB -0.50032611
## GR  0.63616507
## HR -0.23103977
## HU -1.65716237
## IE  0.02869502
## IT -0.46468719
## LT  0.07107211
## LU  0.16907657
## LV  0.33435201
## MT  1.46642982
## NL  0.23828574
## PL -0.04229342
## PT  0.99027518
## RO -0.84670206
## SE -0.75219620
## SI -0.31388036
## SK -0.49254598
## 
## with conditional variances for "country"

# AIC
AIC(glmer_step2)

## [1] 23612.7

# AUC
predicted_probs <- predict(glmer_step2)
auc(scaled_df$trans_docs, predicted_probs)